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Abstract 

The reaction A + B — > B is studied when the reactants diffuse in phase 
space, i.e. their dynamics is described by the Langevin equation. The steady- 
state rate constants are calculated for both the target (static A and mobile 
-B's) and trapping (mobile A and static i?'s) problems when the reaction is 
assumed to occur at the first contact. For Brownian dynamics (i.e., ordinary 
diffusion), the rate constant for both problems is a monotonically decreasing 
function of the friction coefficient 7. For Langevin dynamics, however, we 
find that the steady-state rate constant exhibits a turnover behavior as a 
function of 7 for the trapping problem but not for the target problem. This 
turnover is different from the familiar Kramers' turnover of the rate constant 
for escape from a deep potential well because the reaction considered here is 
an activationless process. 
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I. INTRODUCTION 



Theories of bimolecular chemical reactions in solvents that occur at first contact of the re- 
actants has been traditionally developed assuming that molecules move in diffusive manner^. 
This corresponds to the so-called high friction limit of the more general Langevin dynamics. 
In this paper we develop a theory for Langevin reactant dynamics and analyze the kinetics 
over the entire range of the friction coefficient, from low to high friction. 

Specifically, we study the kinetics of the irreversible reaction A+B — > B in two limiting 
cases: static A and mobile i?'s (target problem) and mobile A and static i?'s (trapping 
problem); the difference between the two problems is illustrated in Fig. 1. Traditional 
theories based on the diffusion equation predict monotonic decrease of the reaction rate 
when the diffusion constant decreases, i.e. when the friction increases. It will be shown that 
in the case of the Langevin dynamics such a behavior of the reaction rate with the friction 
is observed only for the target problem. For the trapping problem the steady state rate 
constant has a turnover behavior as a function of the friction coefficient 7: it first increases 
with 7, reaches a maximum, turns over and then decreases approaching zero as 7 ^ 00. 
This rate constant vanishes at 7 = in one dimension while it remains finite in higher 
dimensions. 

The turnover behavior of the reaction rate is well known in the Kramers' theory of 
activated unimolecular reaction^'!. This theory considers the escape of a Langevin particle 
from a deep potential well. To escape the particle has to overcome a high (compared to keT) 
potential barrier. The escape rate vanishes as 7 — and as 7 00 and has a maximum in 
between. At low friction the rate is limited by energy exchange with the environment and 
vanishes at 7 = since the particle cannot gain the energy necessary to overcome the barrier. 
As 7 becomes finite, the particle can exchange energy with the environment. This leads to 
the increase of the reaction rate. This is the so-called low friction or energy diffusion regime. 
In the opposite limiting case when 7 00 (the high friction regime), particle's motion is 
purely diffusive. As a consequence, the escape rate is proportional to the diffusion constant, 
D = kBT/7 so that the reaction rate vanishes when D goes to zero. 

In the trapping problem a particle does not need to overcome an energy barrier in order 
to react. In fact, even when the friction is zero every particle with finite initial velocity will 
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eventually react. Of course, particles that move faster will be trapped more quickly and 
are the first to be removed from the system. When the friction increases from zero, these 
particles can be replaced by initially slower moving particles that gain energy from the heat 
bath. Thus the trapping rate at first increases as the friction increases. As the friction is 
further increased, the motion of the particle eventually becomes diffusive and in this limit 
the trapping rate is proportional to the diffusion coefficient D which vanishes as the friction 
tends to infinity. 

To analyze the kinetics when motion of the reactants is governed by the Langevin equa- 
tion one has to solve the Klein-Kramers equation in the phase space with appropriate bound- 
ary conditions. This is an extremely complicated task. Even the simple problem of survival 
of a free particle moving on a line in the presence of a single trap has an extremely compli- 
cated solution in the phase space. This solution was found by Marshall and Watsoni forty 
years later after the problem was posed by Wang and Uhlenbecki in 1945. A discussion of 
the source of this complexity and a list of key references can be found in Refs. [ . 

Since we cannot solve the problem in phase space exactly we will find approximate 
solutions in the limiting cases of low and high friction. These solutions are then used to 
construct an interpolation formula that covers the entire range of the friction coefficient. 
The outline of the paper is as follows. In section II we analyze the survival of a Langevin 
particle moving on a line between two static traps. Some of results derived in this section 
has been reported earlieJll without derivation. An approximate solution found in Section II 
is used in Section III to calculate the steady-state rate constant for the trapping problem 
in one dimension over the entire range of the friction coefficient. Since this approach is not 
generalizable to higher dimensions, in Section IV we develop another approximate method 
that can be applied in a space of arbitrary dimension and use it in Section V to find the 
rate constants for the target and trapping problems as a function of the friction coefficient 
in one and three dimensions. Section VI contains a brief summary. 

II. SURVIVAL ON AN INTERVAL TERMINATED BY TRAPPING POINTS 

In this section we analyze the survival of a Langevin particle moving on a line between 
two traps. We calculate the effective rate constant characterizing the kinetics of trapping. 
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By averaging this rate with respect to the length of the interval, we will obtain the steady- 
state rate constant for the trapping problem in one dimension over the entire range of the 
friction. 

Let X and x = v denote the position and velocity of a free particle of unit mass whose 
dynamics is governed by the Langevin equation, 



X + 'jx = R{t) 



(2.1) 



where 7 is the friction coefficient and R{t) is a Gaussian random force of zero mean with 
correlation function given by the fluctuation-dissipation relation, {R(t)R(t')) = 2'y6(t — t'). 
Here and below the thermal energy is set kgT — 1. 

Equivalently, the joint probability density, P{x,v,t), of finding the particle at the phase 
space point (x, v) at time t is described by the Klein-Kramers equation: 



dP{x, V, t) 
dt 



d 



d 



d 



P(a;, f , t) . 



(2.2) 



dx dv dv 

To study the escape of the Brownian particle from the interval [0, 1], we require that P{x, v, t) 
satisfies the absorbing boundary conditions: 



(2.3) 



P{x, V > 0,t) = at X = , 
P{x,v < 0,t) = at X = 1 . 

Assuming that the system has initially a uniform spatial distribution and the Maxwell 
velocity distribution, the survival probability that describes the fate of the particle in the 
interval is given by, 



S{t) 



dx 
^-00 



dv P{x, V, t) 



(2.4) 



where P(x, f , t) is the solution of Eq. (|2.2| ) with initial condition P{x, w , t = 0) = Pcq(a;) feq{v), 
with Peq{x) = 1 and 



feqiv) 



(2.5) 



To characterize the kinetics of the escape from the interval we will use the effective rate 
constant, k, defined as the reciprocal of the mean lifetime of the particle on the interval. 



1 

k 



S{t) dt . 



(2.6) 



When 7 = 0, the particle moves with a constant velocity. The survival probability of a 
particle starting at x with a velocity v is, 



S{t\x, v) = Ri-v) H [x + vt] + Riv) }i[x-l-vt] 



(2.7) 



where H(z) is the Heaviside step function defined as H(z) = for z < and }l{z) = 1 for 
z > 0. Averaging S{t\x, v) over the uniform distribution in x and the equilibrium distribution 



m V gives: 



Sit) 



dx 

"'-oo 



S{t\x,v) /eq(f) dv = erf 
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(2.8) 



This expression shows that in the ballistic limit, S{t) has a power law tail 

1 
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oo 



(2.9) 



From this and the definition in Eq.( |2.6|) it follows that k = when 7 = 0. 

On the other hand, when 7 = 00, the particle does not move and S{t) = 1. As a 
consequence, A; = in this limiting case also. Since k vanishes in the 7^0 and 7 — 00 
limits, this rate constant shows a turnover behavior as a function of 7. 

Since the solution of Eq.( p.2|) with boundary conditions (|2.3| ) is unknown, we adopt the 
following strategy for determining k for arbitrary values of 7. We first derive approximate 
expressions for k in high and low friction regimes, and then use the Visscher-Mel'nikov- 
Meshkovii (VMM) interpolation formula to obtain an expression for the rate constant that 
covers the entire range of 7. 



A. High Friction Regime 

In this section, we analyze the the high friction limit. First, we reduce the Klein-Kramers 
equation ( |2.2| ) with absorbing boundary conditions ( f^.3D to the ordinary diffusion equation 
with radiation boundary conditions. Next, we calculate the effective rate constant by solving 
this equation. 

To proceed, we define the reduced probability density, 

/+00 
P{x,v,t)dv, (2.10) 
-00 
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and the flux, 



r+oo 

j{x,t) = I vP{x,v,t)dv. (2.11) 

-oo 



Integration of Eq.( p.2D over tlie velocity leads to the continuity equation: 

dpjx, t) ^ _ dj{x,t) 
dt dx 



(2.12) 



in which we have used P{x, ±00, t) = 0. Multiplying Eq. ( p.2|) by v and integrating over the 
velocity from —00 to +00, we have: 

dj{x,t) 2 dP{x,v,t) 

— — 7j(a;,t) — / V dv . (2-13) 



dt ' j-00 dx 

This equation can easily be solved to give: 

j{x,t)=e-^'j{x,0) - /We-^(*-*') ^^^^'^'^ dv. (2.14) 

JO J-00 ox 

Note that since the initial distribution of the velocity is Maxwellian, we have j{x, 0) = 0. In 

the 7^1 limit, P{x, v, t') does not vary on times of order of I/7. Therefore, the integration 
over t' can be carried out easily to give: 

./ N 1 /■+°° 9 dP(x,v,t) , , 

j{x,t) = — 1^^^^^^- (2.15 

7 J-00 ox 

Plugging this into Eq.( |2.12| ) leads to: 

ap(.,t)_l r„.3^£i^,„, (2,16) 



dt 7 7-00 dx"^ 

In the 7^1 limit, P{x,v,t) can be written as an expansion in powers of I/7. The zero 
order term of this expansion is, 

Pix,v,t) = f,^iv)pix,t). (2.17) 

Using this in Eq. ( |2.16|) one obtains the ordinary one- dimensional diffusion equation, 

dpix,t) ^ ^ d^p{x,t) . ^2 18) 

dt dx^ ' 7 

To derive boundary conditions supplementing this equation we use the two expressions 
for the flux given in Eqs.( p^.ll]) and ( ^.15| ) to write: 



vP{x,v,t)dv+ / vP{x,v,t)dv 







fO ndP(x,v,t) 
7 J-oo ox 



dv — — 

7 Jo 



^2 dP{x,v,t) 
dx 



dv 



(2.19) 



where the integral over the velocity has been broken into two parts. At the 
boundaries, according to Eq. (p.3|) , P(0,f > 0,t) = P(l,f < 0,t) = 0. This 



suggests that <9P(x, f > 0, t)/9x|^_Q is much larger than (9P(a;, f < 0, t)/(9x|^^Q, and 
dP{x,v < 0,t)/dx\^^^ is much larger than dP{x,v > 0,t)/dx\^^^, when |t>| is not too small. 
Since contribution of small velocities are weakened by the multiplying factor f ^, we neglect 
these small terms. Retaining only major contributions leads to: 



j(o,t) 
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dx 

1 i-o ^dP{x,v,t) 
7 J-oo dx 



dv 



x=0 



x=l 



V P(0, V, t) dv , 

V P(l, V, t) dv . 



(2.20a) 
(2.20b) 



Substituting the approximation for P{x,v,t) given in Eq.( |2.17| ) into these relations, we 
obtain the radiation boundary conditions: 

dp{x, t) 

x=0 

= Kp{l,t), (2.21b) 



D 



dx 
dp{x, t) 



Kp{Q,t) 



(2.21a) 



dx 



x=l 



where k is defined as, 

K={\V\)- 



V\ /cq(f ) dv 



,7r 



2. 1/2 



0.7979 . 



(2.22) 
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Using K given by Eg. ( p. 221) , one can find the corresponding Milne's length, I = D / n = 
1.253/7, which has to be compared with the exact value, / = —(,{1/2) D = 1.46/7 (where 
({z) is the zeta functioni), obtained from a very sophisticated derivation!. The two Milne's 
lengths are in a reasonably good agreement, considering the simplicity of our derivation. 

The survival probability 5'high(^) (where the subscript "high" designates the high friction 
regime) is : 

Shigh{t) = I p{x,t) dx 



(2.23) 

JO 

where p{x,t) is the solution of Eq.( ^.18D with the initial condition p{x,t = 0) = 1 and 
boundary conditions ( |2.21aD and ( p.21b| ). The mean lifetime, k^^^^, of a particle diffusing 



between partially absorbing boundaries and initially uniformly distributed on the interval 
is0: 

^ 7^ + 1^ , (2.24) 



^high ^eq 

where k^q is the rate constant in the reaction controlled limit (i.e., 7 — >■ 0, D ^ oo limit or 
K — i> 0) given by: 

keq = t^PeqiO) + KPcq{l) = 2k . (2.25) 



The second term fc^ in Eg. ( p. 24)) is the rate constant in the diffusion controlled limit (i.e., 



7 —>■ oo, D limit or k oo) which is the inverse of the mean first passage time of the 
particle to the boundaries at x = and x = 1, averaged over the uniform distribution of 
initial positions: 



Plugging this and the expression in Eq.( p.25D into Eq.( |2.24D , we find that the effective rate 
constant is: 

1 1 1 1 7 , , 

+ -— = — + -7- . (2.27) 



fchigh 2k 12D 2k 12 



B. Low Friction Regime 

In this section we analyze the low friction regime where 7 is small and, hence, the particle 
velocity is a slowly varying variable. First, we reduce the Klein-Kramers equation ( p.2|) with 



the absorbing boundary conditions (p.3|) to a diffusion equation along the velocity coordinate 
in the presence of a sink term describing the escape of the particle from the interval. Next, 
we calculate the effective rate constant by solving this equation. 

When 7 = 0, the survival probability of a particle starting at x with initial velocity v is 
given by Eq. { \i-7\j , and the lifetime of such a particle is given by: 

rsaix, v)dt = ^ YLi-v) + iin^ Uiv) . (2.28) 
Jo \v\ V 

Averaging this time with respect to uniform distribution of x gives the mean lifetime, 1 / K{v), 



1 roo 2 

^ , dx S(t\x,v)dt = —- . (2.29) 
K[v) Jo Jo 2\v\ 

K{v) can be regarded as the effective rate constant for the escape of a particle with initial 
velocity v. 

Let f{v, t) be the probability density of the velocity. For finite 7 and in the absence of any 
reaction leading to the depletion of the probability density, f{v,t) satisfies the differential 
equation, 

This equation describes the relaxation of f{v,t) to the Maxwell distribution, feq{v), in 
Eq.( |2.5| ). Because of the escape of particles from the interval, the probability density of the 
velocity decreases with time. The simplest way of taking this decrease into account is to 
incorporate a sink term into Eq.( p.30D as, 

- (2.31) 

where K{v) is chosen so that in the limit 7 — > 0, the mean lifetime of the particle given in 
Eq.( p.29|) is recovered, i.e., K{v) = 2\v\. 

The survival probability S'iow(^) is given by 

roo 

SioAt) = / fiv,t)dv, (2.32) 



where f{v,t) is the solution of Eq.( |2.31|) with the initial condition f{v,t = 0) = fcq{v) 



When 7 = 0, the survival probability obtained by solving Eq.( p.31 ) is. 



/cq(t^)e-^(")*rfi;= / ^=e-2H*rft; = e'*%rfc(y2t) . (2.33) 
00 J —00 \/'2tx 

This function turns out to have exactly the same long time behavior like the exact survival 

probability given in Eq.( |2.9| ). As shown in Fig. 3, the two survival probabilities essentially 

coincide except for short times. 

Integration of Eq.( p.31| ) with respect to v gives the time derivative of Siowif) at t = as: 



dt 



K{v) f,^{v) dv = 2{\v\) =2k = fceq . (2.34) 

t=o 

One can check that dS\Q^/dt = dSiiigi^/dt at t = 0. Another interesting limit of Eq.( p.3lD 



is when 7 — 00. In this case, the velocity distribution instantaneously relaxes to feq{v) 
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so that the probabihty density can be written as, f{v,t) = feq{v) Sio^it), and the survival 
probabihty decays exponentially, Sio^it) = e~'''"^^. From this it follows that the effective rate 
constant is equal to keq in the 7 — oo limit. 



1. Effective Rate Constant 

Let t{v) be the lifetime of the particle initially with velocity v in the presence of the sink 
term, K{v). The effective rate constant in the low friction regime can be expressed in terms 
of r(t>) as 

r+oo 

Kol= r{v)Uq{v)dv. (2.35) 



The lifetime r(f ) satisfies the adjoint equation, 

7e^^/^-^e-^V2^ _2|.;|r(.;) = -l. (2.36) 
dv dv 

The boundary conditions that complement this equation are: r(f ) vanishes as |f | ^ 00 and 
driv) / dv\^^Q = 0, since t{v) is a symmetric function of v, i.e. t(— t>) = t{v). 

When 7 is small, Eq. (|2.36| ) can be solved approximatively by matching the solutions 
found for small and large v. For large velocities, |f | > v*, where v* is a constant to be 
determined, the second term in Eq.( |2.36|) is more important than the first one (since 7 is 
small) and so r(t') ^ l/2|f |. At small v, the second term, proportional to v, can be neglected 
and the exponentials can be replaced by unity. Thus, Eq. (|2.36|) takes the form 



The solution satisfying the condition dr{v) / dv\^^Q = is 



r{v) = r(0) - - . (2.38) 
27 

The time r(0) and the velocity v* are determined from the condition of continuity of t{v) 
and dT{v)/dv at w = v* . This leads to 

l/2v* + (t;*2 - /27 ; Q < \v\ < v\ 
l/2|i;| ; \v\>v*, 
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where v* = (7/2)^/^. 

Averaging t{v) with respect to /eq(f ), one finds that the leading term of A:iow(7) as 7 — > 

is: 

n -1 



k 



low 



3V27r 







In 




m 








2 







h 



(2.40) 



This result is exact and can be derived directly from Eq.( |2.2| ) with appropriate boundary 
conditions. Since /ciow = /ceq as 7 ^ 00, the following heuristic formula can be used to 
interpolate fciow over the entire range of 7: 

-1 

■ , (2.41) 



/ciow = keq <l + — In 
on 



1 + ^ 

7 



in which A is the only unknown constant. It is found that Eq.( |2.41|) accurately fits exact 
results (obtained by solving Eq. (|2.36| ) numerically) for A = 1.45. 



C. Turnover 



So far we have derived expressions for kyow and fchigh as functions of 7. Now, we use the 
VMMii interpolation formula to obtain an analytic expression for k that covers the entire 
range of 7. This leads to 



k{j) _ A;iow(7) ^high(7) 



eq 



1.2 

"'eq 



— In ( 1 + - 

Sir I 7 , 



1 + 



7 



Vl8 



IT 



(2.42) 



In order to test the theory, we have performed Langevin dynamics simulations to compute 
the exact survival probability, S{t), and the effective rate constant, k, as described in Ap- 
pendix The results are reported in Fig. ^ where the closed circles represent simulation 
data, solid and dashed lines give the dependences predicted by Eq.(|2.42|) with A = 1.45 and 



A = 6, respectively, while dot-dashed and long-dashed lines display fchigh (Eq. ( p.27| )) and 
fciow (Eq.( |2.4l|) ) with A = 1.45, respectively. The dotted line designates the purely diffusive 
part of /chigh (see Eq.( p.271 )). 

The rate constant shows a turnover behavior as a function of 7. In the low friction 
regime, it goes to zero like k ~ l/ln(7~^), in contrast to the Kramers' rate constant which 
is proportional to 7. The theory, i.e. Eq.( p.42D with A = 1.45, describes simulation results 
reasonably well. However, Eq. (|2.42|) with A = 6 is in excellent agreement with the simulation 



11 



data as shown by the dashed hne through the data in Fig. 2. Equation ( |2.42|) with A = 6 can 



be regarded as exact for practical purposes. The rate constant k increases as 7 gets larger, 
attains its maximum value (smaller than fceq) close to 7 ~ 1, turns over and then decreases 
like ~ 7~^ as 7 — ^> 00. In this friction regime, the theory is in excellent agreement with 
simulation results. Comparison between the dot-dashed and dotted lines (i.e., diffusion with 
perfectly absorbing boundary conditions) shows how results are improved by using radiation 
boundary conditions. 

In Fig. 3, S{t) is plotted versus t for various values of 7. When 7 = 0, the survival 
probability is strongly non-exponential and goes as at long times (see Eq.( |2.9| )). For 
comparison we have also plotted S'iow(^) given in Eq.( p.33|) . As one can see, 5'iow(^) is slightly 



below S{t) at short times (since the slopes at t = are 2k, and k, respectively), but they 
become indistinguishable from each other at longer times. 

For 7 = 00, S'(t) = 1 since the particle cannot reach the absorbing boundaries. When 
7 = 1, the decay of S{t) is close to exponential but with a slope smaller than k^q = 2k 
represented by e"^"^* (dot-dashed line). As 7 becomes smaller than unity, S{t) becomes 
more and more non-exponential and decays very slowly. When 7 becomes greater than 
unity, S(t) is almost an exponential with the slope /chigh given in Eq.( p.27| ). 

To summarize, the main result of this section is the expression given in Eq. (|2.42| ) for the 
effective rate constant for a particle moving between two absorbing boundaries. It describes 
the turnover behavior of the rate constant as a function of 7. We will use this result in 
the next section to obtain the turnover of the steady-state rate constant in the trapping 
problem. 

III. STEADY-STATE RATE CONSTANT FOR THE TRAPPING PROBLEM IN 

ONE DIMENSION 

Consider an ensemble of traps (absorbing points) uniformly distributed on an infinite 
line. If c is the concentration of traps (i.e., the number of traps per unit of length) and L 
denotes the distance between neighboring traps, the distribution of L is given by, ce~'^^. A 
Langevin particle is injected at any point of the line with equal probability. The probability 
density that a particle is injected on the interval of length L is c^Le~'^^. Let s(t\L) be 
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the survival probability of the particle moving within the interval L averaged over uniform 
distribution of the initial position within the interval and over the equilibrium distribution 
of the initial velocity. Averaging s{t\L) with respect to L, gives the survival probability of 
the particle in the presence of static traps at concentration c: 



POO 

S{t)=c^ Le-'^ s{t\L)dL. 
Jo 



(3.1) 



Before analyzing the turnover of the steady-state rate constant as a function of 7, we briefly 
consider the asymptotic long time behavior of S{t), which is determined by particles located 
in large intervals free from traps. Such a particle changes its velocity many times before 
being trapped. As a consequence, the particle motion is essentially diffusive and its averaged 
survival probability approaches zero as e~^^^^ at very long timeJiil. 

The steady-state rate constant, fcgs, is expressed in terms of the survival probability as. 



rOO /"OO /"OO 

:/ S{t)dt = c^ dt Le-"^ s{t\L)dL . 
Jo Jo Jo 



(3.2) 



Changing variables so as to express s{t\L) in terms of the survival probability on a unit 
interval and carrying out the integration over the time, we have: 



k.Je) 



' z e 



■ dz 



8 = ^ 

c 



(3.3) 



/o k{ez) 

where k{ez) is the effective rate constant given in Eq.( p.42| ). Using this expression we find: 



K 



1=^ hi[ — 

'2-K 37r \e 



97rV27r \e 



(3.4) 



in which the functions hi{y) and /i2(2/) are defined as, 

hr{y) = 27E + 21n(y) - y + [y'' - 2y + 2] eJ^Ei(y) , 

h2{y) = 67E + 6 \n{y) - Ay + y^ - ^ - 3y^ + 6y - g] e^ E,{y) 



(3.5a) 
(3.5b) 



where 7e = 0.577216 is the Euler's constant and, Ei{y) = J^z~^ dz, is the exponential 
integral functionll. 

Equation ( p. 41 ) predicts a turnover behavior of the steady-state rate constant as a function 
of e. When e ^ 0, kss approaches zero as 



k.Je) 371 



K 



B = A exp 



37r - 3 + 27E 



(3.6) 
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In the opposite limit when e — > oo, the steady-state rate constant tends to zero hke, 



K 



(3.7) 



In Fig. 4 the ^-dependence of A;ss(£^) given in Eq. (|3.4|) is shown for A = Q (solid curve) and 
A = 1.45 (dashed curve). The former may be considered as exact since it is obtained with 
A = 6, which provides an almost perfect fit to simulation data for the unit interval (see Fig. 
2). The dashed curve shows the dependence obtained on the basis of the theory developed 



in Section [II B| . Although both curves have the same qualitative behavior there is clearly a 
room for improvement. In what follows we develop another approach which allows one to 
calculate kssie) in a space of arbitrary dimension. We will use the dependence shown by 
the solid curve in Fig. 4 in order to test the approach in one dimension. It will be shown 
that the new theory predicts the dependence of kss{s) which is in better agreement with the 
exact result (solid curve in Fig. 4) than does the theory developed in Section |II Ij| (dashed 
curve in Fig. 4). 



IV. NEW APPROACH 



We now introduce a new and simpler approach to calculate the steady-state rate con- 
stant for the trapping problem which has the great advantage of being generalizable to any 
dimensions. Moreover, this approximate approach is expected to improve as the dimension 
of the space increases. As we will see below, even in one dimension, the new approach is 
already in better agreement with the simulations than does the previous one (compare Figs. 
4 and 5). 

The basic idea of this new approach is to apply the VMM interpolation formula directly 
to the steady-state rate constants found in the low and high friction regimes. In the low 
friction regime the rate constant will be obtained by using a many-particle generalization 
of the approach presented in Section [II B[ In the high friction regime we will assume that 
the steady-state rate constant for the trapping problem is equal to the rate constant for the 
target problem. This approximation is quite good in one dimension and is virtually exact 
in three dimensions. 
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A. Target Problem 



The survival probability of the target particle (i.e., the A in the reaction A + B — > B 
with Da = 0), is given by 

5targ(t) = exp |-c^* A;targ(t') dt'^ , (4.1) 

where c is the concentration of mobile traps (i.e., the 5's) and fctarg(^) is the time-dependent 
rate coefficient. The steady-state rate constant i 

K.l, = c St.,,{t)dt. (4.2) 

The exact time-dependent rate coefficient /ctarg(^) is obtained by solving the Klein-Kramers 
equation with appropriate boundary conditions. This is an extremely complicated problem^ 
and here we will approximate /ctarg(^) by the CoUins-Kimball rate constant, kcj<i{t), which is 
obtained by solving the diffusion equation with the radiation boundary condition at contact, 
in which the intrinsic rate constant is chosen so that A;ck(0) is exact. 

First, we note that kta.rg(t = 0) can be found exactly for any friction. The product 
cfctarg(O) is just the collision frequency, therefore 



{\v\) = J2/71 ] d=l, 
nRmvl) = VS^R^ ; d = 3, 

where R is the contact radius. In addition, when the friction vanishes, 7 = 0, fctarg(^) = 
^targ(O) and thus the rate constant is given by Eq.(O) for all times. As the friction tends to 
infinity, the rate constant is given by the Smoluchowski theory. Now, we consider the Collins- 
Kimball rate coefficient, /cck(^), obtained for the diffusion constant £) = I/7 and the intrinsic 
rate constant k appearing in the boundary conditions equal to fctarg(O) given in Eq. (|4.3| ). 
This Collins-Kimball rate coefficient has the following properties: i) kcj<i{0) = A;targ(0) at 
any friction, ii) when the friction vanishes, we have A;ck(^) = ^targ(O) for all times, and in) 
at high friction kcK{t) coincides with the Smoluchowski rate constant except for very short 
times where it is equal to fctarg(O). In the view of these properties we approximate fctarg(^) 
by kcK(t). As a consequence, we have Stargit) = ScK{t) and 

KJ.g = c ScK{t)dt. (4.4) 
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B. Trapping Problem in the Low Friction Regime 



When the friction vanishes, the particle moves with a constant velocity. Its survival 
probability is equal to the probability that the region visited by the particle in time t is free 
from traps. Since the traps are uniformly distributed in space, the probability that a region 
is free from traps is given by the Poisson distribution. The survival probability of a particle 
moving with a velocity v is 

{exr)(—c\v\t) : d = 1 , 

(4.5) 
exp{-CTcR^\v\t)) ; d = 3. 

Assuming equilibrium distribution of the initial velocity, the probability density of the ve- 
locity at time t is 

f{v,t) = S{t\v)f,^{v). (4.6) 

Integrating this with respect to velocity, we find the survival probability in the ballistic 
regime 

5baii(t) = / f{v, t)dv= I S{t\v) /eq(t;) dv . (4.7) 



This result is exact. In one dimension it leads to 

^baiK^) = y-j exp |-cft - y I c/t; = exp f — j erfc f j . (4.8) 

This expression can also be obtained using S{t) in Eq.( p.8|) in the relation for S{t) given in 
Eq.( |3.lD . In three dimensions 5baii(^) is 

1/2 



/2\V2 .00 r ^2^ 

Shaiiit) = i-j exp I -CTiR^vt - ^ \ dv 



^1 + '—^ exp erfc r J - - ; r = Kci = cRH , (4.9) 

where we have used the relation in Eq. (|4.3 ) for k in three dimensions, i.e. k, = y/Sir R'^. 

Equations ( [4.5| ) and ( [4.6| ) show that particles with large velocities react more rapidly and 
therefore the higher the velocity the faster the probability density decays. When the friction 
is non-zero, due to the interaction with environment, the probability density, depleted by 
the reaction, relaxes to the Maxwell distribution. At small friction f{v,t) satisfies 
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^i^ = [C{v)-k{v)]f{v,t), (4.10) 
where C{v) describes the relaxation of the velocity distribution in d-dimensions 

(4.n) 

and the sink term k(v) describes the reaction. In one and three dimensions the sink term is 

elf I ; d = 1 , 
k{v) = { (4.12) 

TiR^ c\v\ ; d = 3 . 

When 7 = 0, C{v) vanishes, and f{v,t), obtained from Eq. (|4.10|) with initial condition 
/(f , 0) = feq{v), coincides with Eq.(^]B|). 

The steady-state rate constant can be obtained by averaging the mean lifetime t{v) of 
the particle with the initial velocity v, with respect to the equilibrium distribution of the 
initial velocity: 



oo /-oo 



k^ol = I ^'''feqi^)r{v)dv ■ v'-'f,^{v)dv = l. (4.13) 

The lifetime t{v) can be determined by solving the backward equation 

'C^{v) - k{v)] t{v) = -1 , (4.14) 

where C^iv) is the adjoint operator 

C\v) = 4- e^^/2 -^f e-"^/2 4- ■ (4-15) 
^ ' f^-i dv dv ^ ^ 

In summary, to find the steady-state rate constant for the trapping problem one has to 
calculate the rate constants fciow and /chigh corresponding to low and high friction regimes 
and then use the VMM interpolation formula. The rate constant kiow is given by Eq.( [4.13D 
where r(f ) is obtained by solving Eq. ([4.14|) . The rate kugh is assumed to be approximatively 



equal to the steady-state rate constant for the target problem fctarg given in Eq. ([4.4{). In the 
following section we use this approach to calculate the steady-state rate constant over the 
entire range of friction in one and three dimensions. 



17 



V. STEADY-STATE RATE CONSTANT 



A. One-Dimensional Case 



We begin with the calculation of fc^arg based on Eg. ( [4. 41 ) using the survival probability 
<ScK{t) = exp < -c 



'iGDt AD 



'kH\ 



TX 



K 



exp I ^ 1 erfc 



Id 



(5.1) 



where D = l/'-f and k is given in Eq. 

fctargfe) I VStT I VStT 



K 



exp 



As a result, we obtain 



^ + e^erfc(v/^) -1 | rf^j ; 



6=- 
C 



(5.2) 



The ratio kts,rg{e)/K, as a function of e is represented by the dash-dotted line in Fig. 5. One 
can see that fctarg(£^) / monotonically decreases with e from unity at e = to zero as £ ^ oo. 
For sufficiently large e the steady-state rate constant approaches the one predicted by the 
Smoluchowski theory and tends to zero like 1/e. 

To calculate fciow by Eq. (|4.13|) we have to determine r(f) solving Eq.( [4.14| ) which in one 
dimension is 



7e 



,,2 Ir, d 



dv 



dT{v) 
dv 



C\V\T[V) 



-1 . 



(5.3) 



By changing variables this equation can be reduced to Eq. ( |2.36| ). This allows us to use the 
solution in Eq.( p.4lD and to eventually obtain 



i^low ) 



1 + ^lnfl + A) 
Stt V 2eJ 



(5.4) 



where A = 1.45. The dashed line in Fig. 5 represents k\o^{e)/ k as a function of e. The ratio 
monotonically increases with e from zero in the ballistic regime {e = 0) to unity as £ — oo. 
The rate constant /ciow(^) vanishes as l/ln(l/e) as e — 0. 

We now use the VMM interpolation formula to get the steady-state rate constant for the 
trapping problem for the entire range of friction (i.e. of the parameter e = '~f/c) 

^trap(^) _ ^low(^) ^high(^) _ ^low(^) ^targ(^) 



(5.5) 



Using the expressions in Eqs. (|5l^) and (|5.4|) we find 
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^trap(£) 



1 + — In 



1 + 



1.45 
27 



57r 



exp 



/Stt 



^ + erfc (v/^) 



(5.6) 



This function is represented by the long dashed hne in Fig. 5. It should be compared with 
the exact solution found in Section |T| which is represented by the solid line in Fig. 5. The 
approximate result in Eq.( ^.6|) is overall closer to the exact dependence than the solution 
found in Section |T| (see Fig. 4 and related discussion). However, Eq. (|5.6|) is no longer exact 
in the diffusive limit {e —>■ oo). This is because the steady state rate constant of the trapping 
and target problems differEl. However, this difference is negligible in three dimensions. 



B. Three-Dimensional Case 



Like in the one-dimensional case, we use the CoUins-Kimball survival probability in 
Eq.( [4.4|) in order to calculate ktarg- In three dimensions 



<ScK{r) = exp 



l + e 



T + 



a 



+ e"^ erfc (v^) - 1 j 



r = CKt 



(5.7) 



where k is given by Eq.(^]3|) and 



a 



'l+e) 



(5.8) 



where = 47r_R^c/3 is the volume fraction of traps. Substituting i5ck(t) into Eq. ([4.4|) we 



obtain 



^targ(^; ' 



■ POO 

/ ScK{T)dT 

Jo 



(5.9) 



In three dimensions the steady-state rate constant is a function of two variables e and 
For small As </>, Eq.( ^.9|) reduces approximatively to 



^targ (^) ' 



l + e 



\ (1 + 



(5.10) 



Next we find the mean lifetime t{v) by solving Eq. ([4.14|) which in three dimensions has 
the form 
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7 ^.72^ 
2 dv 



V 



driv) 
dv 



— nR cvt(v) = —1 



(5.11) 



Introducing an effective diffusion constant X = •y/cn = 2TTe/3(j), we can rewrite Eq.( |5.11| ) as 



2 dv 



V 



driv) 
dv 



VT[V] 



1 



(5.12) 



Tliis equation can be solved for small and large A. When \ ^ oo (instantaneous relaxation 
of the velocity distribution) t{v) = I/ck for all v. In the opposite limiting case of small 
A the first term in the right-hand side of Eq. (|5.12| ) can be neglected except for very small 
values of v. However, when v is small the exponentials in the first term in Eq.( p.l2| ) can be 
replaced by unity. Thus at small A, t{v) satisfies 

2 driv) 



A d 

v"^ dv 



dv 



— vt{v) = . 

° CK 



(5.13) 



A solution of this equation that is finite at f = and vanishes as f — oo is 



t{v) 



Ai (7rVV2i/2^i/3 



Ai(0) 



CUV 



(5.14) 



where Ai{z) is the Airy function^. 

Once r(t>) is known one can find fcio^ by Eq.( ^.13D . Using r(f) in Eq.( |5.14| ) we find that 
as A — 0, fciow is given by 

A;iow(A) vr , 3^/6 



[r(2/3)]' A'/^^ . 



(5.15) 



K 4 47ri/3 

This expression shows that in three dimensions fciow remains finite in the ballistic regime 
(A = 0), in contrast to the one-dimensional case where fciow vanishes as friction goes to zero. 
Since A = 27r£:/30 = 7/cK, Eq.( p.l5| ) shows that the rate constant increases as 7^/^ for small 
7- 

As A 00, t{v) = 1/ck, for all v and the ratio k\o^{e)/ becomes unity. To interpolate 
between the two limiting cases of small and large A we use the formula 



fclow(A) 



f + aA2/3 



3^/6 [r(2/3)]' 



3.643 



(5.16) 



K l + aA2/3 ' 7rV3(4-7r) 

To test this formula, we calculate A^iow/'^ by numerically solving Eq.( p.l2D and then carrying 
out the integration with respect to v in Eq.( [4.13| ). The result is reported in Fig. 6 which 
shows good agreement between the dependence predicted by Eq. (|5.16|) (solid line) and the 



20 



numerical data (closed circles). The dashed line shows the dependence predicted by Eg. ( ^.15|) 
for small e. 

To cover the entire range of friction we again use the VMM interpolation formula 

/Ctrap(£, 0) _ fclow(A) fctarg(£, 0) 



(5.17) 



where k\o^{\)/K and fctarg(^, 0)//^ are given in Eqs.( |5.16D and (^]9|), respectively. When 
<^ 1 , the second term in the right hand side of Eq. ( |5.17| ) becomes independent of and is 
simply given by the first term in Eq. (|5.10| ). Thus, in the -C 1 limit, Eq. (|5.17|) reduces to 



/^trap(£,0) _ l+h{e/4>f^ 



« (1 + e) [l + 6(5/0)'/' 

where 



(5.18) 



, /27r\2/3 3i/6(47r)i/3[r(2/3)]' , , 

U^^'^^ ^ 5.964. (5.19) 

Figure 7 displays the dependence of ktr:a.p{e , (p) / k (solid lines) as a function of e for two 
values of 0. The dashed and dot-dashed lines represent fciow/ ^ and fctarg/ ^, respectively, for 
the same 0. The rate constant exhibits a turnover behavior as a function of e, however, the 
turnover in three dimensions depends on the volume fraction and is much less pronounced 
than in one dimension. 



VI. SUMMARY 

In this paper we have generalized the standard theory of diffusion controlled reactions to 
the case when the reactants diffuse in both coordinate and velocity space (i.e., they undergo 
Langevin rather than Brownian dynamics). We have developed an approximate theory of 
the steady state rate constant for both the target and trapping problems. This theory was 
tested against accurate results obtained from simulations for the trapping problem in one 
dimension and is expected to work even better in higher dimensions. 

The key finding was that the steady state rate constant for the trapping (but not for 
the target) problem exhibits a turnover behavior as a function of the friction coefficient for 
Langevin dynamics. For Brownian particles the rate constants for both problems decrease 
monotonically as the friction increases. The physical explanation of this turnover is that for 
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the trapping problem in the baUistic regime, a particle with near zero velocity can survive 
for a very long time. In one dimension, where the most probable velocity is zero, the mean 
lifetime is actually infinite and thus the rate constant is zero in this limit. Increasing the 
friction increases the rate because particles with initial velocities close to zero can be speeded 
up by random forces. In three dimensions, where the most probable velocity is finite, the 
rate constant is also finite in the ballistic limit and hence the turnover is less pronounced. 
Since we have assumed that the particles react on first contact, there is no energy barrier to 
reaction. Thus the friction dependence of the trapping steady state rate constant represents 
a simple physical example of turnover in activationless rate processes. 



22 



APPENDIX A: THE SIMULATION PROCEDURE 



Simulations were performed using the discretized version of Eq. 

Xn+l = Xn + e~'^'^ {Xn - Xn-l) + Xn , (A.l) 

where A is the time step and the Gaussian random noise X„ (related to i?„ by, X„ = 
[A(l — e''^^)/'-]'] Rn) is defined by the moments, 

= and (X„X„,) = (^^^ (1 - e-^^)^ 5^^, . (A.2) 

For each trajectory, the initial position xq is generated from the uniform distribution between 
and 1 and the initial velocity Vq from the centered Gaussian distribution of unit standard 
deviation. From this, the first position xi at the first step is then calculated as: 

xi = Xq + Vq^ , (A. 3) 

and the next positions x„ are generated according to the algorithm in Eq.( [A.l| ). To simulate 
the absorbing boundary conditions, each trajectory starting at Xq (0 < Xq < 1) at time t = 
is terminated at time tj = nA when either the condition x„ < or x„ > 1 is satisfied for the 
first time. The first passage time and the survival probability Si{t) (defined as Si{t) = 1 
for all t < ti and Si{t) = otherwise) for this given trajectory are recorded. The survival 
probability, S{t), and the effective rate constant, k, (i.e., the inverse of the mean lifetime) 
are then obtained by averaging over a large number of trajectories: 

Sit) = ^ES^it) and k-' = -J2u. (A.4) 

i=i i=i 

For all simulations reported in this paper we used the time step A = 10'^ for 10"^ < 7 < 10^, 
and N = 10^ trajectories were used to perform the averages. 
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FIGURES 
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FIG. 1. Sketch of the reaction A + B — ^ B between a single A with many B's where the 
circles show the contact radius for the reaction. In the target problem (panel I) a single A is 
fixed in space and surrounded by mobile B's while in the trapping problem (panel II) the S's are 
immobile and the A is mobile. In this latter case the static traps B's are distributed according to 
the Poisson distribution and overlaps may occur. 

1 I I I — I I I I — I I I I — I I I I — I I I I — I I I — I I I I nil 
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FIG. 2. Reduced effective rate constant, k/keq, as a function of the friction coefficient 7. The 
data (closed circles) are obtained from Langevin dynamics simulations, the dot-dashed line corre- 
sponds to /chigh given in Eq.( |2.27] ), the long-dashed line to /ciow given in Eq.(|2Hl), the solid and 
dashed lines through the data represent the VMM interpolation in Eq.( |2.42| ) with A = 1.45 and 
^ = 6, respectively. The dotted line is kd/keq, where kd, given in Eq.( p.26D , represents the purely 
diffusive part (i.e., k 00) of fchigh- 




10.0 



FIG. 3. Survival probability S{t) as a function of time for various values of 7. Solid lines 
correspond to the Langevin dynamics simulations results. For 7 = 0, the solid line represents the 
exact result in Eq.( |2.8| ) and the dashed line Sio^{t) given in Eq.( p.33 ). The dot-dashed line is. 



-2Kt 
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FIG. 4. Reduced steady-state rate constant, A:ss(e)/K, as a function of the dimensionless friction 
coefficient e = "f/c. Solid line (Eq.( |3.4D with ^ = 6) is the essentially exact result while the dashed 
line (Eq. (|3.4| ) with A = 1.45) is obtained using the approximate theory. 
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FIG. 5. One-dimensional reduced steady-state rate constant, A:trap(£)//«, for the trapping prob- 
lem as a function of the dimensionless friction coefficient e = 7/c. The solid line represents the 
exact result (Eq.(^) with A = 6), the dot-dashed line corresponds to A;targ(e)/'« given in Eq.( |5.2[ ) 
for the target problem, the dashed line to k\ovi{^)/i^ given in Eq.(5.4) with A = 1.45, and the 
long-dashed line represents the VMM interpolation in Eq.(^.6|). 
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FIG. 6. Three-dimensional reduced rate constant, k\o^{\) / in the low friction regime as a 
function of the parameter A = 27re/3(/). The closed circles are obtained from numerical solution of 



Eq. (|5.12[) , the solid line corresponds to expression in Eq. (|5.16D and the dashed line to the small-A 
behavior in Eq.( |5.15| ). 
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FIG. 7. Three-dimensional reduced steady-state rate constant, ktra,p{^,<P)/f^, for the trapping 
problem as a function of the dimensionless friction coefficient e = "yR/^/^jr for two values of the 
volume fraction of traps (p. The dot-dashed line corresponds to fctargC^j </')/'^ given in Eq. (|5.9| ) for 



the target problem, the dashed line to ki^w given in Eq.( 5.1q ), and the solid line represents the 
VMM interpolation in Eq.( |5.17|) . 
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